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We present a novel formulation to calculate transport through disordered superconductors con- 
nected between two metallic leads. An exact expression for the current is derived, and is applied 
to a superconducting sample described by the negative-L'' Hubbard model. A Monte Carlo algo- 
rithm that includes thermal phase and amplitude fluctuations of the superconducting order pa- 
rameter is employed, and a new efficient algorithm is described. This improved routine allows 
access to relatively large systems, which we demonstrate by applying it to several cases, including 
superconductor-normal interfaces and Josephson junctions. The effects of decoherence and dephas- 
ing are shown to be included in the formulation, which allows the unambiguous characterization 
of the Kosterlitz-Thouless transition in two-dimensional systems and the calculation of the finite 
resistance due to vortex excitations in quasi one-dimensional systems. Effects of magnetic fields can 
be eeisily included in the formalism, and are demonstrated for the Little-Parks effect in supercon- 
ducting cylinders. Moreover, the formalism enables us to map the local super and normal currents, 
and the accompanying electrical potentials, which we use to pinpoint and visualize the emergence 
of resistance across the superconductor-insulator transition. 
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I. INTRODUCTION 

Chief amongst the remarkable effects observed in su- 
perconductors is their eponymous perfect conductivity. 
Within BCS theory [Ij, where superconductivity arises 
due to pairing between electrons, the effects of temper- 
ature r, magnetic field B, and disorder are well under- 
stood: as the pairing amplitude is suppressed by these 
physical parameters, the system becomes normal, and 
attains a finite resistance. For low-dimensional systems, 
on the other hand, it has been long understood that 
phase fluctuations of the pairing amplitude play a ma- 
jor role in the loss of perfect conductance In two- 
dimensional systems, for example, it has been demon- 
strated [7 that as the temperature increases there is a 
critical temperature Tkt where vortices and anti- vortices 
unbind and proliferate through the system, leading to 
the loss of global phase coherence and superconductivity, 
even though the pairing amplitude remain finite. Indi- 
cations of such a Berezinsky-Kosterlitz-Thouless (BKT) 
transition have been observed in Josephson-junction ar- 
rays [1], in superconducting (SC) thin films [5 , and pos- 
sibly in high-Tc cuprates [^. 

In recent years there has been a reinvigoration of re- 
search into low-dimensional superconductors. This has 
been motivated by intriguing experimental observations 
of electronic transport through disordered SC thin films, 
such as a huge magnetoresistance peak [Tj and a "super- 
insulator" phase [ff, and by the technological progress 
in producing two-dimensional superconductors in the in- 
terface between two oxides [3] and in making ultra-thin 
cuprate superconductors T^. Many of these observa- 
tions are not yet satisfactory explained, chiefly because 
there is no theory that can calculate the current, even 
numerically, through a disordered superconductor, based 



on a microscopic model. 

The calculation of the resistance within the BCS 
picture, usually based on the Bogoliubov-de Gennes 
(BdG) mean-field approach, is straightforward. Blon- 
der, Tinkham, and Klapwijk (BTK) [Tl] studied the 
reflectance and transmission at a metal-superconductor 
junction, and an analogous study was performed 
at superconductor-metal-superconductor junctions |12) . 
Similar approaches [13] utilized the Buttiker-Landauer 
picture |141 115j for non-interacting Cooper pairs to study 
scattering through a SC region. (A difficulty with the 
direct application of the BdG formalism is the non- 
conservation of charge, which can be overcome by study- 
ing a normal ring containing a SC segment jl6j.) An 
alternative approach near to the BCS critical tempera- 
ture is to use a scaling assumption for the conductiv- 
ity [17j . The current through diffusive normal metal- 
superconductor structures has also been calculated using 
a Keldysh scattering matrix theory (TS]. All these ap- 
proaches neglect phase fluctuations so cannot be used 
to study two-dimensional superconductors that exhibit a 
BKT-like transition at low temperatures. 

The resistance of low-dimensional superconductors can 
also be calculated using phenomenological models. The 
conductivity of uniform systems can be probed analyt- 
ically by studying phase slips across the sample within 
the Ginzburg-Landau approach [19l [20] . Thermally ex- 
cited phase slips explain both non-linear conductivity 
and vortex creep induced resistance [21j , whilst quantum 
activated phase slips can drive SC wires insulating |22j . 
However, phenomenological calculations are neither un- 
derpinned by a microscopic model nor include Coulomb 
repulsion or disorder except for the introduction of a phe- 
nomenological normal state resistance. 

Here we develop a new formalism to calculate the cur- 
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FIG. 1: (Color online) A schematic of the experimental setup 
within the negative-?/ Hubbard model. The left and right- 
hand metallic leads are shown in blue, from which electrons 
can tunnel through the barriers shown by the gray links into 
the central SC region which is shown in red. 



rent through a superconductor taking into account phase 
fluctuations in the presence of disorder, finite T and B, 
and Coulomb repulsion. The approach we detail here is 
based on the Landauer-Buttiker scheme [HI [15] , where 
one attaches metallic leads to the sample, and then calcu- 
late its conductance. The lead-superconductor tunneling 
barriers ensure that the conductance of the system is al- 
ways finite, even in the SC phase. A previous attempt 
using the quantum Monte Carlo approach to calculate 
current in disordered systems employed the fluctuation- 
dissipation theorem via the current-current correlation 
function |23j . 

The Landauer formula [HI [TS] is a widely adopted 
method to calculate the current through a mesoscopic 
sample that contains non-interacting particles. Meir and 
Wingreen [21] have generalized the formula to produce an 
exact expression for the current through any interacting 
region attached to non-interacting leads, which has been 
successfully applied to a wide range of systems. Follow- 
ing this approach, we partition the system into the three 
parts shown in Fig. [T] the left-hand lead, the central 
interacting region, here a superconductor, and the right- 
hand lead. In the leads the natural particle basis set are 
electrons, and in the sample the natural basis set are Bo- 
goliubons. To circumvent this mismatch of particle basis 
sets we reformulate the Meir- Wingreen formula in a Bo- 
goliubon basis set to derive an exact expression for the 
current flow through a possibly SC region, attached to 
two metallic leads. 

Having derived a general, exact formula, the SC re- 
gion is then modeled by a generalized negative-?/ Hub- 
bard model based on the lattice shown in Fig. [T] In- 
troducing two local auxiliary flelds (which reduce to the 
local density and gap at zero temperature), we decouple 
the interacting fermions. While the conductance formula 
is exact, in order to evaluate correlation functions we 



neglect quantum fluctuations, and integrate numerically 
over the thermal fluctuations of the auxiliary parameters 
[551 [2S] using a Monte Carlo method. A signiflcant ad- 
vantage of the formalism is that it allows us to construct 
current and potential maps of the system. These allow 
us to diagnose the microscopic features that increase the 
resistance of the sample. This paper details the new pro- 
cedure and presents a number of applications of the for- 
malism for simple systems, where one can compare with 
existing theories. 

The paper is organized as follows: in Sec. |ll] we flrst 
derive an exact expression for the current through a SC 
region. Using this expression, in Sec. HI we describe how 



the current can be calculated numerically and outline im- 
provements to the auxiliary field approach that allows us 
to study large systems. Having developed the new for- 
mula for the current and accompanying computational 
tool, it is vital to carefully test it against a series of known 
results. Therefore, in Sec. |IV A] we study superconductor- 
normal interfaces in clean systems, and compare with the 



BTK transmission formulae, while in Sec. IV B we study 



the temperature dependent current in a Josephson junc- 
tion. In Sec. IIV CI we describe how effects of decoher- 
ence and dephasing are manifested in the formalism. We 
investigate the temperature dependence of resistance in 
Sec. |IVD| in which we uncover the temperature depen- 
dence of the resistance and the nonlinear I — V behavior 
that characterizes the BKT transition in two dimensions 
and vortex excitations in quasi one-dimensional systems. 
We then, in Sec. |IVE[ apply an external magnetic field 
to probe the Little-Parks effect. Finally, we demonstrate 
how to construct current and potential maps for the sys- 
tem, and use them to study the microscopic behavior 
at the superconductor-insulator transition in Sec. |IVF| 
The details of the analytical derivation and the numerical 
procedure are described in the appendices. 



II. ANALYTICAL DERIVATION 
A. Current Formula 

To calculate the current for interacting particles we 
start with the general formula for the current [23] 
through an interacting region, connected between two 
non-interacting leads 

/ de[Tr{(/L(6)rL-/R(e)r«) {Gl - GD] 



(1) 



Here /^(e) = [exp(/3(e - /i^)) + 1]"^ with x e {L,R} is 
the Fermi distribution of the left (L) and right-hand (R) 
leads that are held at chemical potentials and reduced 
temperature j3 = \/k^T (where Ub is the Boltzmann 
constant). The imposed potential difference eV = A/i = 
A'L — Mr between the leads drives the current J through 
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the system. The integral is over all electronic energies e. 
^ij = 27rl]aex'^«('^)^'»*^<y channels a in lead x, and 
Ya.i is the tunneling matrix element from channel a in 
the the lead to site i in the sample. Finally, Glj^, Qfjai 
and Qf-^ are the electronic retarded, advanced, and lesser 
Green functions (in the site basis) for electrons of spin a 
in the sample calculated in the presence of the leads. 

Eqn. ([I]) is exact, and captures, via the electronic 
Green function Q, all the processes that can transfer 
an electron through the system. When the intermedi- 
ate regime has SC correlations, some of these processes 
involve Andreev scattering - absorption of an electron 
pair by the condensate and a propagation of the remain- 
ing hole. To expose these processes, it is convenient to 
transform from the electron basis set {cl^,Cia) with site 
index i into the Bogoliubov basis set {'^na-Tlna), using 
the Bogoliubov-de Gennes relations Cia- = Ui{n)^na~ 
av* {n)jn_„ (at present Ui and Vi are arbitrary, except for 
the unitarity condition, but later on they will be deter- 
mined by the actual Hamiltonian that will be used for the 
SC region). The Green functions transform from the elec- 
tron basis into the energy basis set of Green functions 
{G^,G^} and the family of anomalous Green functions 
H>(m,n) = -i(7m-a7L)> H<_(to,7i) = , 
H>(m,n) = -i(7m-o- Trier), and H<(m,n) = ii^i^crlma) 
according to 

= u, (G> - G<) u* + V, (G>, - G< J V* 

- av* (H> - H<) u* - an, (H>^ - H< J vj , (2) 



and 



^<(i,j)=u*G<u;-v,G>,v*+au*H>v*-av,H5,u,. 



(3) 



Solving for the Green functions across the system in the 
presence of the leads (Appendix A), leads to the final, 
exact result 



J 



-y 
h ^ 



de[/L(e)-/R(e)]x 



Tr 



(p - R iG'^n.'^.H'' -f (r^. . - w .)Gi'''r"'^H^'' 

V uv ' y\s)^(7 ' V* ' 'cr ' V u* v* ' v* u* /^cr ' uv ' 'cr 

ari. H^(r-^.-r-4. )G^+ar^,. Ht-(r,->^- r:/)Gt'- 
^(r^.. + r,V)(H^r-^.H; + Ht-r^/Ht;)! , (4) 



where rjy(m,n) = 27r^ 



now in the transformed basis set. This is written in 
a form describing transmission from the left-hand side 
to the right-hand side of the sample. We will show in 
Sec. |IV A| that it therefore exposes the rise of resistance 
due to the suppression of correlations between the left 
and right-hand sides of the superconductor. 



We note that deep in the SC regime where the SC gap 
obeys Y, and in the case where the leads inject elec- 
trons within the gap such that eV < 2 A, we can make a 
perturbative expansion in small tunneling Y . This yields 
the simple expression for the current 



J 



--T 



Tr 



de[/L(e)-/R(e)]x 

[(r^., + rj.jG:(r,-^- 



(5) 



Here G are the Green functions calculated in the absence 
of the leads. This equation has direct Y"^ dependence 
on the tunneling matrix element, with neglected higher 
order contributions of order ^ (y/A)^, as it describes 
Cooper pairs tunneling through the contact barrier. We 
shall show later that this contribution is precisely what 
is predicted for the current |22i according to the BTK 
formula [11], and notably, as the leads inject electrons 
only into the gap, there is no normal current, but only 
Andreev processes allow the flow of current. The pertur- 
bative form Eqn. ([s]) offers two important computational 
advantages. Firstly, it is considerably less resource inten- 
sive to calculate as the Green functions are diagonal so 
it does not demand summations over separate variables. 
Secondly, it does not require the expensive matrix inver- 



sion embodied in Eqn. ( A2 ) to find the general equation 
for the current. Due to its usefulness we also note that 
an analogous expression can be derived for the normal 
current when injecting electrons outside of the gap 



J 



-T 



deTr[/L(e)(r^., + r„^.J 

-/R(e)(r?.u + r^*v)]^^G:., (6) 



where 3 stands for the imaginary part. Since this term 
represents the normal current, it has a direct Y^ depen- 
dence on the tunneling matrix element. Though they 
offer a considerable computational advantage, these per- 
turbative formulae cannot be used on the border of the 
superconductor-insulator transition where the supercon- 
ductor gap breaks down and A < y. Therefore, un- 
less specified, we use the full expression for the current, 
Eqn. Q, in our numerical calculations. 



B. Current and voltage maps 

Eqn. (|6|, with a coefficient of Y^ , describes the nor- 
mal current that enters and leaves the system as single 
electrons, whereas Eqn. ([sj with a coefficient of Y'^ corre- 
sponds to a tunneling supercurrent. However, the normal 
and supercurrent can interchange inside the sample. In 
order to understand the microscopies behind phenomena 
in the disordered superconductor it is vital that we can 
probe the spatial distribution of the current as it switches 
in nature through the sample. Therefore, here we extend 
our formalism to map out the flow of current within the 
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sample. To calculate the current distribution map we 
use the general expression for the current crossing a sin- 
gle bond [551 US] from site i to j 

J^^ = ^iT.J^ [t^,S<{J,^) - t,,g<{^,J)] . (7) 

(7 

Transforming again into the diagonalized basis, the local 
current is 



2ev-^ fde 
^ h 2n 

I <-/ / *v u 'u V 



A'i -AH 

u* u 'v* V 



G<- 



uu* vv* 



H< + cr fA'^'^ - A'-' l H 



a f ' 



(8) 



where t\]l^{m, n) — tij\ii{m)w j{n) . As before, the normal 
and anomalous Green functions are calculated in 
Eqn. (A3) in the presence of the leads. Moreover we 



note t 



;hat the current comes in two flavors, the contri- 
bution to the current from the normal Green function 
is associated with the normal current and that from 
the anomalous Green function H< gives the Cooper pair 



current. In Sec. IV F we verify that this intersite current 
yields the correct net conservation of charge. 

To provide an additional probe into the nature of the 
superconductor-insulator transition we extend the for- 
malism to map the local chemical potentials across the 
sample. This should reveal any weak links and the lo- 
cation of the sources of resistance in a sample. To de- 
termine the local effective potential at a specific site we 
add a weak link from that site to a third lead (a "tip"). 
The tunneling current from the tip into the sample is then 
calculated, and the chemical potential of the tip adjusted 
until that current flow is zero. This chemical potential 
thus corresponds to the effective local chemical potential 
in that site. To calculate the current flow into the tip we 
first evaluate the full Green functions G in the sample in 
the presence of voltage drop between the left and right 
reservoirs Eqn. ( A2 ) but without the tip. We then use the 
perturbative formula for the current, Eqn. ([5]), but with 
one lead representing the left/right hand leads, and the 
other the perturbative tip. This process is repeated for 
each site in the sample (due to the perturbative nature 
of the tip, this calculation can be done simultaneously 
for all sites). In Sec. IV F we demonstrate how maps of 



the potential can expose weak links in the sample and 
help diagnose the microscopic mechanisms that give rise 
to resistance. 



III. MODEL AND NUMERICAL PROCEDURE 

In the previous section we have developed an exact for- 
mula for the current through an arbitrary intermediate 
region, which may include SC correlations. We now use a 
specific model to describe this SC region - the negative-?/ 
Hubbard model, a lattice model that includes on-site at- 
traction, and may include disorder, orbital and Zeeman 



magnetic fields, and even long-range repulsive interac- 
tion (which we will not deal with in this paper). The 
Hamiltonian is 

-^Hubbard — ^ ^ ^iac\^Cicr ^ ^ t^icj^^cj^^Ci^Ci^ 

-E ' (9) 

where eia- is the on-site energy, tij the hopping element 
between adjacent sites i and j, and Ui is the onsite two- 
particle attraction, taken to be uniform, i-independent, 
in this paper. An orbital magnetic field can be incorpo- 
rated into the phases of the hopping elements tij , while 
a Zeeman field splits the spin-dependent on-site energies 
Ci^. In this paper we will only deal with orbital fields. 
To account for disorder, will be drawn from a Gaus- 
sian distribution with characteristic width W . The inter- 
site spacing is a. Unlike, for example, the disordered 
XY model, the negative-C/ Hubbard model can lead to 
a BCS transition, a BKT transition, or to a percolation 
transition, and thus this choice is general enough not to 
limit a priori the underlying physical processes. Impor- 
tantly, the model includes the fermionic degrees of free- 
dom which may be relevant to some of the experimental 
observations. 

Calculation of correlation functions, for example the 
Green functions that enter the current formula, require 
thermal averages. To perform the thermal average we 
need to decouple the quartic interaction term so we em- 
ploy the exact Hubbard-Stratonovich transformation 



VAVAe 



- fi dr -I%^+A,(r)cl,ct^+A,(r)c„c,t 



(10) 



which is basically a Gaussian integration {T>A = 
Ilr.idAi{T), where the product runs over all times and 
all sites). Note that the field Ai{T) is just an integration 
variable that decouples the two-body term in the super- 
conducting channel, and should not be confused with 
|[/i|(ci^Ci^). Similarly, one introduces the integration 
fields pirj{T), that couple to the spin density [30] {c\^Cia), 
and leads to an additional term — ^ \Ui\pi-a{T)<^i„Cia 
in the action (which, in the mean-field approximation 
gives rise to the Hartree-Fock contribution) [50,. Decou- 
pling in both the A and p channels not only provides ac- 
cess to both soft degrees of freedom, but also the saddle 
point solution gives the standard mean-field results for 
those fields, and furthermore guarantees that the action 
expanded to Gaussian order corresponds to the random 
phase approximation |31| . 



The Hubbard-Stratonovich transformation ( 10 ) is ex- 
act. Since our main interest lies in thermal effects, for 
example the thermal BKT phase transition, or thermal 
activation of vortices, we now neglect quantum fluctua- 
tions (the r dependence of A). One can now write the 
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partition function for the Hubbard model as 



Tr 



-^-f^Hubbard 



2?(A,p)Trf 



-^«BdG(A,p) 



(11) 

where the latter trace is over all fermionic degrees of free- 
dom. HBdci^, p) is the Bogohubov-de Gennes (BdG) 
Hamiltonian with a given set of A and p, where these 
vectors designate the set of values of these parameters on 
all lattice sites 



(12) 



Given the explicit form of the diagonalizable BdG 
Hamiltonian, we can calculate expectation values and 
correlation functions, 



Tr 



pO 



N 



2?(A, p)e-''^« ^ ( 



o 



(13) 

where the sum is taken over all positive eigenvalues 
(quasi-particle excitations) of the BdG Hamiltonian. 
Here Eq, and \n) are the ground-state energy, excita- 
tion energies and excitation wave functions, respectively, 
for the BdG Hamiltonian, for the specific configuration 
of A and p. It is straightforward to see that in this case, 
the saddle-point approximation of the partition function 
gives rise to the mean-field BdG equations (and then A.; 
indeed corresponds to \Ui\{ci^Ci^)). The calculation of 
the full integral, using the (classical) Monte Carlo ap- 
proach [35], includes also the contributions of thermal 
fluctuations of the amplitude and phase of the order pa- 
rameter. In Appendix [B] we detail how we improve on 
contemporary methods to perform the Monte Carlo cal- 
culation in 0{N^^^ IvP/^) time, where N is the number 
of sites and M the order of a Chebyshev expansion. 



IV. APPLICATIONS 

We have derived a new expression for the current flow 
through a superconductor, and demonstrated how to cal- 
culate the current in mesoscopic systems. Before apply- 
ing it to understand and predict novel phenomena, it is 
important to verify it across a variety of exemplar sys- 
tems, where one can compare against well-established 
theories. As the main novelty of the approach is the 
inclusion of thermal fluctuations, we pay particular at- 
tention to verifying the formalism in two dimensions, 
especially looking for signatures of the BKT transition 
driven by phase fluctuations. In Sec. |IV A] we probe the 
current through a clean superconductor, and check that 
we can recover the BTK results for the contact resis- 
tance. A key effect in such systems is Josephson tunnel- 
ing, so in Sec. |IV B we study the temperature dependence 
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FIG. 2: (Color online) The lower graph (b) shows the varia- 
tion of conductance a with interaction strength. The normal 
current is shown in red and the supercurrent in blue. The 
numerical results are shown with error bars, and BTK theory 
by the solid line. The upper plot (a) shows the variation in 
the order parameter with interaction strength; numerics are 
shown with error bars and the mean-field theoretical predic- 
tion by the solid black line. The horizontal green line denotes 
the chemical potential difference across the sample. The gray 
shading denotes the range of interactions where < A < eV , 
so that both normal and Andreev processes contribute to the 
current. 



of the resistance of a single Josephson junction. De- 
phasing (by temperature averaging) and decoherence (by 



electron-electron interactions) are studied in Sec. IV C 



In Sec. IV D we examine the temperature dependence of 
the resistance on the two sides of the BKT transition and 
compare to analytical results. In Sec. |IVE| we introduce 
flnite magnetic fleld (flux) and probe the Little-Parks ef- 
fect in the presence of disorder. Finally, in Sec. IV F 
we demonstrate how plotting maps of the current and 
potential across the system can illuminate the micro- 
scopic processes at the superconductor-insulator tran- 
sition. Throughout we use an attractive interaction of 
U — 1.6t to describe the SC region. To avoid the Van 
Hove singularity at half filling (33! we study systems at 
an average 38.7% filling, except for Sec. IV C where we 
focus on wires with a low filling fraction of 20%. In the 
linear response regime we impose a potential difference 
of eV = 0.02t. Systems were typically two-dimensional, 
so a single lattice site thick, 12 lattice sites wide, and 48 
lattice sites long. 



A. Clean systems 

At low temperatures, thermal fluctuations of the pair 
amplitude and phase may be neglected. The conduc- 
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tance, in this limit, through a clean SC region has been 
calculated by Blonder, Tinkham and Klapwijk (BTK) 
[TT] , and is solely due to the contact resistance at the two 
interfaces. By assigning a tunneling strength XjZ to the 
barriers, BTK have shown that if the intermediate sample 
is in the normal state, then the current is purely due to 
electrons tunneling across the barrier, and the transmis- 
sion coefficient is given by l/Z^ Jlj. On the other hand, 
if the sample is in the SC state, then the SC gap, A, in- 
hibits electrons from directly tunneling into it. Instead, 
these electrons Andreev tunnel accompanied by a hole. 
For a large barrier Z ^ 1 the transmission coefficient 
becomes A^/4Z''(A^ — E^) [TT], where E is the electron 
energy. Electrons with an energy outside of the gap can 
either tunnel alone with a corresponding normal trans- 
mission coefficient (E + ^ E'^ - t^^) j {2Z'^ E'^ - A2), or 
Andreev tunnel with accompanying hole, and have a 
transmission coefficient of I^Z^{E'^ — A^). We first 
compare the results of our numerical calculations to these 
BTK formulae, and then demonstrate that for the simple 
case of a single SC site, the BTK results can be derived 
analytically from our current formula. 

To verify that the model recovers the correct behavior 
at the tunneling barrier we focus on the weak coupling 
limit. In this limit, once a Cooper pair tunnels through 
the first barrier, it has an equal probability of continu- 
ing to either the left or the right lead, an consequently 
the current through the double barrier will be half that 
of a single barrier [27 . For a long enough system the 
finite bias and temperature smear any Fabry-Perot type 
interference. To study the effect of the changing order 
parameter A, we focus on a 39% filled system at "zero" 
temperature (without quantum fluctuations) , where A is 
indeed equal to the pair correlation |J7|(ci^Ci^), vary the 
interaction strength J7, and monitor the various compo- 
nents of the tunneling current. For a pristine system with 
= 0, all of the resistance stems from the two tunneling 
barriers, and we verified that the current flow was inde- 
pendent of the length of the SC region. A relatively large 
potential bias of eV = O.lt was applied across the leads. 
This allows us to explore all tunneling processes, either 
for A < eV or A > eV by changing the interaction pa- 
rameter U and as a result A, see Fig.|2][a). Our results for 
the current are depicted in Fig. [2]jb)7 and has Z k, 1.4. 
At J7 = the current is entirely normal. As shown in 
Fig. [2ja), with increasing U the SC gap grows giving rise 
to a resonance in the Andreev current when = A. At 
the same time the normal current falls as fewer electrons 
can be directly injected outside of the SC gap. As the 
interaction strength is increased further, so that the SC 
gap exceeds the chemical potential difference, resonant 
electrons are no longer injected into the divergent den- 
sity of states at the SC gap, and the Andreev current 
falls. In agreement with the BTK calculation, at large 
A the Andreev current adopts its final value, 1/4 of the 
normal ?7 = conductance and no normal current flows. 
The agreement with the BTK prediction verifies that the 
current formula Eqn. (pi) contains the correct tunneling 



behavior. 

We now turn to derive the BTK results from our for- 
malism analytically, which can be done straightforwardly 
in the weak coupling limit y ^ 1 when the SC region 
consists of a single site, and we take the leads as having a 
parabolic dispersion. In the linear response regime where 
a potential V is put across the sample such that injected 
electrons are entirely within the SC gap, we find that the 
normal current is zero and we recover the analytic result 
for the Andreev current, J = e^V t^^ j^hZ'^^t^^ - M^), 
where Z — ^ ^/ttv/Y , /x is the chemical potential, and u 
is the density of states at the Fermi surface. If the sample 
is normal we find that there is no Andreev current, and 
the normal current is J = e^V /2hZ^ . These results are 
what would be expected from the BTK formalism, and 
coupled with the numerical results confirm that the for- 
malism properly treats tunneling between the leads and 
the SC sample. 



B. Josephson junction 

Another simple example that we wish to explore is 
a single Josephson junction, which will be modeled in 
the negative-C/ Hubbard model by an intermediate region 
consisting of two clean superconductors, between which 
we insert a weak link where the nearest-neighbor hop- 
ping element t' is smaU {t' « t), see Fig. [sj Studying 
this system will allow us to probe how a phase difference 
across a barrier can affect the current flow through it. 
We again adopt a 39% filled band with no disorder. We 
first set i' = to disconnect the left and right-hand sides, 
and numerically evaluate the current through the central 
region. This is by no means trivial. The current formula, 
through the anomalous Green function, allows an absorp- 
tion of a pair from the incoming lead into the condensate 
on one side of the barrier, and an emission of another 
pair into the outgoing lead. This current, however, will 
depend on the phase difference between the SC order pa- 
rameter on the two sides of the barrier. For t' = 0, i.e. 
an infinite barrier, the phases of the left and right-hand 
order parameter are uncorrelated, and thus all phase dif- 
ferences are degenerate in energy. Therefore, the current 
vanishes, but only after averaging over all states, which 
is done automatically in our numerical procedure. In the 
other limit, when the hopping matrix elements are the 
same as the hopping through the rest of the supercon- 
ductor, t' ~ t, the phase of the superconductor is locked 
so we see the standard free SC current flow. 

We now model the situation with a moderately sized 
central barrier. This splits the superconductor in two, 
but crucially a Josephson supercurrent flows between the 
two sides, thus allowing the current to flow with no addi- 
tional resistance. The current J(T) = Jj(r) cos(0l — 0r) 
is maintained by the phase difference 0l — 0r between the 
left and right-hand superconductors, and the maximum 
value of the dissipationless current is the critical Joseph- 
son current Jj(T) = (tt] A|/2ei?n) tanh(| A|/2fcBr) [M], 
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FIG. 3: (Color online) Upper: The setup to model the Joseph- 
son junction. Traversing the center of the SC region is a 
Josephson junction (opaque cuboid). The junction is modeled 
by the reduction of the matrix hopping elements to t' = 0.005 
(brown interconnects) compared to f = 1 in the supercon- 
ductor. The two metallic leads are shown in blue, and the 
lead-superconductor tunneling barrier by the gray cuboids. 
Lower: The variation of conductance with temperature for 
the Josephson junction. Results of the numerical computation 
(points) and of the theoretical model (solid lines) are shown 
for a weak (green), intermediate (red), and strong (blue) cou- 
pling between the two superconductors. 



where i?n is the resistance of the central barrier when 
the system is in the normal state. 

In order to study the thermally driven disruption of the 
Josephson current numerically, it is vital that this break- 
down occurs before the BKT transition occurs, which as 
we show in Sec. |IV D[ by itself reduces the current flow 
through the system. We thus use a small hopping ele- 
ment for the tunneling barrier of t' = 0.005i, which has 
a large pn and therefore small Josephson current Jj. In 
Fig. [3] we show the current as a function of tempera- 
ture, in the presence of the weak link. When there is 
no voltage drop across the Josephson junction, the cur- 
rent Jm that flows through it is given by Jm = V/R, 
where R is the contact resistance to the normal leads. 
This current is maintained as long as the critical Joseph- 
son current Jj is larger than Jm- As temperature is in- 
creased thermal fluctuations will weaken the phase lock 
between the two superconducting regions, and the crit- 
ical current is reduced. When Jj is reduced below Jm, 



a finite voltage develops across the Josephson junction. 
This drives the phase difference across the junction to 
increase with time, which in turn leads to an oscillating 
current. This current has a non-zero time- average [T], 
leading to a total resistance Rjj{T) = R/{1 — -s/l — A^), 
where A = Jj / Jm < 1 [I] ■ This time averaged current is 
exactly the quantity calculated in our Monte Carlo pro- 
cedure. Fig. [3] depicts a comparison between this simple 
model and our full numerical calculation, with reasonable 
agreement. 

The critical current can be modified by varying the 
resistance of the central barrier, R^- The intermediate 
case has i?n — 0.075/i/e^, the stronger coupling with 
Rn — 0.06/i/e^ is obtained by lowering the barrier to t' ~ 
O.Oli, and the weaker coupling with R^ = 0.085/i/e^ by 
widening the original barrier {t' = 0.005i) to four lattice 
sites. This wider barrier weakens the coupling between 
the superconductors so the Josephson resistance emerges 
at a lower temperature. A lower barrier strengthens the 
coupling so raises the temperature required for the emer- 
gence of resistance. Both these regimes are consistent 
with our simple model. We also verified that at very 
strong coupling where the temperature required for the 
breakdown of phase coherence becomes of the order of 
the BKT transition temperature our simple model for 
the current fiow no longer captures the full physics of the 
system. This study validates that our formalism can cor- 
rectly model the presence of the Josephson supercurrent 
across the weak link introduced into the superconductor, 
and can therefore be used to model mesoscopic systems 
that contain multiple SC grains. 



C. Decoherence and dephasing 

The issue of decoherence and dephasing plays a sig- 
nificant role in transport at low temperatures. Here we 
define decoherence as the many-body phenomenon that 
leads to the loss of coherence via interactions among the 
electrons or interactions with the environment. On the 
other hand, dephasing can occur in a non-interacting 
system, and emerges from the fact that due to the fi- 
nite temperature, electrons possess a range of energies, 
of the order of ksT. Electrons of different energies ac- 
quire different phases along their respective trajectories, 
and if these phases differ by 27r or more when their en- 
ergy changes by fceT, then interference phenomena will 
average out to zero. 

Decoherence: The effects of decoherence due to 
electron-electron interactions are more profound in one- 
dimensional wires in the normal phase. Since the original 
Hubbard model employed in this calculation (Eqn. ^) 
is an interacting model, one expects decoherence to arise 
naturally from the calculation. However, though the orig- 
inal formula for the current is exact, the approximation 
employed above - Eqn. ( 13 1 - does not include quantum 



fluctuations. This means that it neglects the imaginary 
component of the self energy which corresponds to damp- 
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FIG. 4: (Color online) (a) The relative fall in current in a 
one-dimensional sample due to the introduction of self energy 
at k^T — O.lfi. The black points show the numerical results, 
and the red line highlights the expected theoretical variation 
with length [3^, where Jo is the current that flows when 
impeded solely by the contact resistance, (b) The relative 
change (A J = Jo — J) in current as a function of temperature 
in a one-dimensional normal phase sample. The black points 
show the numerical results, and the red the expected model 
variation. The two green dashed lines show the exp{— jj, / ksT) 
and behavior, (c) The changing current in the presence of 
a SC phase in a two-dimensional sample. The vertical green 
dashed lines show the BKT and normal phase transitions. 

ing due to interactions, and the resulting decoherence. In 
order to test the effects of such decoherence due to many- 
body interactions, we introduce into the normal Green 
function for momentum k (as here we study a wire in the 
normal phase) , by hand, the self energy 

lim y- ^(Cp)[l-^(ep-q)][l-»(ek+q)] .^4. 

5^o27r4-^ w + Cp-^p-q-a+q-i^ ' 

which is the lowest order contribution to the single- 
particle self-energy, and where are the momentum en- 
ergy eigenstates of the Hamiltonian. 



A similar approach [36" has been applied to interacting 
electrons in a continuous one-dimensional system with 
repulsive contact interactions (the second order contri- 
bution to the self-energy does not depend on the sign 
of the interaction). In this case, it has been shown, for 
wires with parabolic dispersion and chemical potential /i, 
that this damping leads to a change in the distribution 
function and reduction in the conductivity by a factor of 
1 - TT'^{k^T/^ifL/l2[L + eexp{^j,/kBT)] where the 
wire length, L, is long enough that the smearing of the 
Fermi surface due to scattering (that occurs over the re- 
laxation length scale £ [32]) outweighs that due to tem- 
perature. For sufficiently long wires L ^ iexp{^/kBT) 
the reduction in conductivity becomes length indepen- 
dent l-TT^{kBT/ll)^/12. 

In order to be able to compare to this theory (which 
relies on the parabolic dispersion), we focus on a system 
with a low filling fraction of 20%, near the bottom of the 
band, and set the disorder to W = O.lt. We employ the 
perturbative expression for the current, Eqn. ([6]), to give 
us access to long wires with L ^ £. In the upper panel of 
Fig.|4]we show the fall in current, as a function of length, 
due to the inclusion of self energy at k^T = 0.07/i, here 
£ « 27a. The overall change of ~ 0.5% is small due to 
the Pauli blocking of scattering processes near the Fermi 
energy. We see reasonable agreement with the model 
over a range of length scales. The middle panel of Fig. |4] 
depicts the change in current with temperature for a sys- 
tem of a fixed length. We highlight the agreement to 
the expected variation in the fall in conductance with 
temperature |36| . At low temperatures {k^T <C /i) the 
damping is severely Pauli blocked so the characteristic 
damping length-scale exceeds the system length and the 
current correction l — TT^{kBT/ii)^exp{—^/kBT)L/12£ is 
exponentially suppressed. As temperature increases the 
Fermi liquid behavior starts to dominate the correc- 
tion to the current. At high, usually unphysical temper- 
atures (k^T ^ /i), numerics see a smaller current shift 
than predicted by theory as the details of the specific 
Hubbard band dispersion versus the parabolic dispersion 
in which the model was developed become important. 

In Fig. |4|^c) we examine the effect of a SC phase on 
decoherence. At low temperature the presence of the SC 
gap suppresses many-body scattering processes. How- 
ever, when temperature is raised above the BKT phase 
transition, scattering events are possible, though have a 
smaller impact on the current than in the normal phase, 
due to the still finite local pair correlations. Above the 
mean-field BCS phase transition the current follows the 
expected parabolic profile as in the normal phase. Thus 
we have demonstrated that while quantum fluctuations, 
as they affect decoherence, can be taken into account in 
our formalism, their effect on the current, for the range 
of parameters studies here, is usually small at < 1%. We 
are thus justified in neglecting them in this study. 

Dephasing: Having observed decoherence in the sam- 
ple we now turn to study dephasing due to thermal av- 
eraging. To verify that our formalism captures this im- 
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FIG. 5: (Color online) The fall in conductance with length 
in a non-interacting one-dimensional wire with disorder W = 
0.2t at two different temperatures. The red trend lines show 
a linear drop off in conductance with length, and green an 
exponential decay. 



portant phenomenon, we study the length dependence 
of the conductance in non-interacting systems. We first 
verified, for a non-interacting clean system, that the net 
macroscopic current increases by 2e^//i for each new con- 
duction channel introduced (not shown), independent of 
length. Setting the amplitude of the disorder toW = 0.2t 
and working at 39% filling, in Fig. [5] we show the fall in 
conductance with length at two different temperatures. 
At T = there is an initial linear fall in conductance over 
length scales smaller than the localization length ^ w 93a, 
and an exponential fall at greater lengths. This is in 
accordance with the expectations of Anderson localiza- 
tion [37 - for length scales below the localization length, 
the conductance changes as a power law of the length, 
while it decays exponentially when the length becomes 
larger than the localization length. At k^T = O.Oli, on 
the other hand, dephasing causes different parts of the 
system to be incoherent with respect to the others, caus- 
ing the conductance to fall linearly with inverse length, 
as expected from a classical system. This observation 
confirms that the formalism naturally incorporates the 
physics of dephasing in disordered systems. 



D. Variation of resistance with temperature 

We have now verified that our formalism captures 
the basic phenomena of contact resistance in Sec. |IV A[ 
Josephson coupling in Sec. |IVB[ and dephasing and de- 
coherence in Sec. |IV C With these key tests complete, 
we are now ideally poised to study further effects within 
the superconductor, starting with the temperature de- 
pendence of the conductivity and its relation to the 



BKT transition. With increasing temperature a two- 
dimensional superconductor undergoes a BKT transi- 
tion [35] characterized by the emergence of vortices across 
the system, leading to the loss of global phase coherence. 
At a higher ("mean-field") temperature, the SC order 
is completely suppressed and the system loses the SC 
correlation even locally. To study how this transition 
is reflected in the current flow we performed numerical 
simulations on a two-dimensional 39% filled SC system 
at several different temperatures. Simulations were per- 
formed for two different levels of disorder, W = 0.1t and 
W — 0.2i, to determine how the transition and current 
flow are modified by the normal-state resistance, and ex- 
trapolated over length to remove the effects of the contact 
resistance (Fig. [6]jb)). 

Even at temperatures below the BKT transition, vor- 
tices can be nucleated from the edge of the sample and 
traverse the system, driven by the Magnus force due to 
the finite current. This produces dissipation at any non- 
zero temperature and current J, according to the non- 
linear potential V cx J^+'^t^t/t ^ Fig.jeja) we see 
that below the BKT temperature the linear resistance, 
that is limy-j-o V/J^ is zero. The plots Fig. |6][d) show 
several simulations that were performed for different im- 
posed potential differences V across the sample, which 
allowed us to extract the index 7 of the conductance rela- 
tion F cx JT. In Fig.[6];c) we show that the conductance 
relation approximately follows the expected theoretical 
behavior with 7 = 1 2T^t/T. 

At temperatures above the BKT transition vortices 
and anti-vortices can easily unbind, though they may be 
partially pinned by disorder. The finite conductance G 
of a sample in this case has been shown by Halperin and 
Nelson 21j to be given by 



G = 0.37G„(e+/ec)' 



(15) 



where Gn is the normal state conductance, is the SC 
coherence length, and ^+ is the SC order correlation 
length, which diverges at Tkt- The critical behavior at 
temperatures near the BKT transition T > Tkt leads to 
the conductance 



G = 0.37Gn6-^ exp[VKrc-TKT)/(r-rKT)] , (16) 

where 6 is a number of order unity. At temperatures 
higher than the (renormalized) mean-field critical tem- 
perature Tc, the conductance is given by the Aslamasov- 
Larkin theory [HI [3S] 



G = 0.37G„(Te - TKT)/(r - Tkt) • 



(17) 



Finally, we can also estimate the crossover between these 
two regimes by noting that the difference between the 
Kostcrlitz-Thouless and mean-field transition tempera- 
tures critical regime is given by |21j 



(18) 



The difference between these two temperatures therefore 
widens with falling normal state conductance. 
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FIG. 6: (Color online) (a) The variation of resistance with temperature for two different values of disorder calculated numerically 
(points). The red solid line shows the theoretical low temperature behavior, and the blue solid line the theoretical high 
temperature behavior. The dashed vertical green lines show the BKT Tkt and mean-field Tc temperatures, (b) The numerical 
results (black points) and the deduced linear length dependence (red line) of the resistance for three different points (i) T < Tkt, 
(ii) Tkt < T < Tc, and (iii) T > Tc- (c) The current nonlinearity index 7 in V oc for two disorder levels, as a function 
of temperature. The black points are the numerical results and the red lines are from theory. The vertical green dashed line 
highlights the BKT transition temperature, (d) Examples of the numerical measurements of dimensionless current against 
voltage leading to the values of 7 shown in the left-hand plot. The best-fit lines employed are shown in red. 



In Fig. [6ja) we depict the variation of resistance with 
temperature above the BKT transition, showing the two 
types of dependence on temperature as is expected by 
theory. We also note that the rising disorder increases 
the normal state resistance CTn, and also broadens the dif- 
ference between the Kosterlitz-Thouless and mean-field 
transition temperatures, which agrees with Eqn. ( 18 ) 
within 20%. The high temperature Aslamasov-Larkin 
expression for the conductance persists well above the the 
mean-field critical temperature, where the SC state has 
been totally suppressed. Finally, when the temperature 
is of the same order as the bandwidth, fc^T ^ t, the resis- 
tance in Fig.|6j[a) increases superlinearly as the Fermi dis- 
tribution becomes smeared across the whole band struc- 



ture. 

Having studied the nonlinear J — V characteristic 
in two dimensions we now turn to look at the one 
dimensional system. Here thermal fiuctuations can 
drive the formation of phase slips at any temperature 
and so this system has the J — V characteristic V = 
Joi?sinh(J/Jo) [111201110], where Jo = 4:ekBT/h and 
R = (/i/4e2) X {Hn/kBT)exp{-AF/kBT) is the resis- 
tance with attempt frequency hfl fa 3. It and energy bar- 
rier AF sa 3.7i. In Fig. [T^a) we show the consistency of 
the numerical model both for the nonlinear J — V char- 
acteristic, and in Fig. [7]jb) the variation with tempera- 
ture. The strong accord between analytics and numerics 
in both one and two dimensions gives us confidence that 
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(a) Nonlinear J — V 




(b) Resistance 




FIG. 7; (Color online) (a) The nonlinear J ~V characteristic 
of a one-dimensional superconductor at fixed temperature. 
The numerical points are shown with black error bars and the 
Langer-Ambegaokar-McCumber-Halperin model [191 [20| by 
the red line, (b) The variation of resistance with temperature 
for fixed bias. The numerical points are shown with black 
error bars and the Langer-Ambegaokar-McCumber-Halperin 
model by the red line. 

the formalism can be applied to study and explore less 
well understood mcsoscopic superconducting systems. 



E. Little-Parks effect 

Varying an applied magnetic field has long been an 
important experimental probe of the properties of a su- 
perconductor. It is therefore imperative to verify that 
the current formula developed here, coupled with the 
Hubbard model for the superconductor, is able to ac- 
curately model the effects of an applied magnetic field. 
In the Hubbard model the effects of the magnetic field 
are incorporated, via the Peierls substitution, into the 
phases of the hopping elements, Uj — )■ ii^e^'^"^'^/'*" where 
00 — hcje is the quantum flux, and the phases 4>ij are 
defined such that their integral over a closed trajectory is 
equal to the magnetic flux threading the surface spanned 
by the trajectory. 

In order to check whether this procedure captures the 
effect of an orbital magnetic field, we apply it to a hol- 
low cylindrical superconductor, of radius r, such as that 
shown in Fig. [Sj threaded by magnetic flux. As demon- 
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FIG. 8: (Color online) Upper: A schematic of the cylindri- 
cal wire within the negative-?/ Hubbard model. The left and 
right-hand metallic leads are shown in blue, from which elec- 
trons can tunnel through the gray toroids into the central 
SC region which is shown in red. The magnetic flux thread- 
ing the cylinder is shown in green. Lower: The variation of 
current with longitudinal magnetic flux at T = Tc. The com- 
putational points with error bars are shown in black, and the 
Little-Parks model best fit is shown by the red dashed line. 



strated by Little and Parks [41 , the flux suppresses su- 
perconductivity and the transition temperature falls pe- 
riodically with the flux. This is often probed by measur- 
ing the falling conductance of the cylinder near to the 
transition temperature [HI 02] . 

We apply our formalism to the cylindrical thin-walled 
superconductor shown in Fig. [s] at 39% flUing and no 
disorder, and apply an external magnetic flux (p along 
the axis of the cylinder. The additional phase shift to 
the hopping matrix elements around the cylinder cir- 
cumference causes the energy of electrons in the cylin- 
der of radius r to increase with trapped flux (p as 
h^{n -f 20/(/)o)^/2r7ir^, where the integer n is chosen to 
minimize the energy. This results in a periodic parabolic 
variation of the electron energy with flux and thus a 
parabolic periodic oscillation in the SC transition tem- 
perature ATc = h^{n + 2(l)/(l)oy /l&mr'^ gT]. Therefore, 
for a cylinder held just below its superconducting tran- 
sition temperature, with increasing flux the supercon- 
ducting state is disrupted periodically and the resistance 
varies with flux, as a series of parabolas, with minima in 
the conductance at every half flux quantum 4> = n4ni/2. 
This has indeed been observed experimentally [41j . 

In Fig. [8] we take a cylinder held near to its SC tran- 
sition temperature and numerically evaluate the conduc- 
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FIG. 9: (Color online) (a) The average fractional error in 
conservation of current X^i^i \^Ji\/ JiN on each site against 
the fraction of total states K/N included in the calculation of 
the current, (b) The changing conductance (black line) with 
width I/norm of Central normal region, the right axis shows the 
normal fraction Jnorm/ Jtotai of the total current (blue line) 
flowing through the central region. 
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tance as a function of the magnetic flux. The reasonable 
agreement with theory demonstrates that the formalism 
correctly picks up the effects of an applied magnetic field. 
The deviation from the parabolic predictions of mean- 
field theory at every half flux quantum is due to thermal 
fluctuations, and will elaborated upon in a later publica- 
tion. 



FIG. 10: (Color online) The upper panel shows the potential 
difference V{x) across the sample with total potential drop 
V. The lower panel shows current maps for short (a) and 
long barriers (b) respectively. Supercurrent is shown by cyan 
darts and normal current by violet pointers, arrow length cor- 
responds to current magnitude and orientation to the direc- 
tion of current flow. Color density corresponds to the order 
parameter |A|, which has peak value Aq. 



F. Current distribution maps 

One important feature of our formalism is the new ca- 
pability to map out the flow of both super and normal 
currents within a sample and the changes in chemical po- 
tential which drive that flow. Since we can now study the 
current flow around impurities in the sample and expose 
weak links with large potential drop, we should be able 
to probe phenomena in the disordered superconductor 
with unprecedented detail and trace their cause back to 
a microscopic mechanism. While applications of this for- 
malism to the outstanding problems in this field will be 
described in future publications, in this section we aim 
to demonstrate the usefulness of the current and poten- 
tial maps, first by further studying the Josephson junc- 
tion with a superconductor containing a central normal 
region, and secondly by studying the superconductor- 
insulator transition in disordered systems. However, we 
will first verify our current mapping formalism by exam- 



ining the site-by-site current conservation in a 39% filled 
system with no disorder. As the only sources and sinks of 
current are the two metallic leads, a consistent calcula- 
tion should obey charge conservation for all of the inner 
sites of the sample. In Fig.jojja) we show the average frac- 
tional error in conservation of current X^iLi \^Ji\/ on 
each site as we vary the number of states K included in 
the calculation out of a possible N states, as prescribed 
in the penultimate paragraph of App. |Bj We see that if 
only 5% of states are included there is a 20% average 
leakage of the current. However, if we include 50% of the 
states in the calculation of the current there is a leakage 
of only ~ 2%. Throughout the remainder of this section 
we include 40% of the states in the calculation to yield 
an average error of approximately 3%. 

Having verified the conservation of current, we demon- 
strate what can be learned from the current maps by first 
studying a modified Josephson setup consisting of two 
clean 39% filled SC regions with a central normal region 
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(b) SC side of transition, T ^ O.UTc 
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(c) At the superconductor-insulator transition, T ~ Tc 
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FIG. 11; (Color online) (a) Shows the fall in conductance 
across the superconductor-insulator transition. Current maps 
on tuning temperature from (b) a superconductor at T ~ 
0.14Tc through to (d) an insulator at T ^ 2.3Tc. At T ^ 
the superconductor-insulator transition takes place. Super- 
current is shown by cyan darts and normal current by violet 
pointers, arrow length corresponds to current magnitude and 
orientation to the direction of current flow. Color density cor- 
responds to the order parameter |A|. Lines of equal chemical 
potential are shown in white. In the current map (b) three 
points of interest are labeled: (1) the normal state, (2) the 
superconductor state, and (3) Josephson tunneling. 



that has U = 0. We can then monitor the current flow 
through the system to see it change from SC to normal in 
character as the intermediate normal region is widened 
in Fig. [ojb). For a narrow U — central region the two 
SC regions are phase locked and predominantly a Joseph- 
son current flows (lower panel in Fig. [lO{^a)). Due to the 
strong proximity effect, the system is entirely SC with 
no reduction in conductance. The electrical potential is 
dropped on the two contact barriers, and remain constant 
through the superconductor (upper panel in Fig. [lOj[a)). 
(For the present case of two equal contact barriers the 
potential in the SC is equal to the average of the chemi- 
cal potential of the two leads). On the other hand, when 
the central U = region is wide, £norm ^ 4a, the two 
SC regions are too weakly coupled for a Josephson cur- 
rent to flow, and instead a normal current flows between 
the two SC regions (lower panel in Fig. [Tojjb)). This, in 
turn, introduces a new resistor into the sample and the 
conductance drops accordingly. Now the potential drop 
is mostly across the Josephson junction (upper panel in 
Fig. [lO|^b)) - the left-hand superconductor adopts, ap- 
proximately, the potential of the left-hand lead and the 
right-hand superconductor that of the right-hand lead. 
This situation is analogous to current flowing between 
SC grains in a disordered sample, and can reveal whether 
they are coherently coupled, when a supercurrent flows 
between the grains, or decoupled, when a normal current 
flows. Such analysis could be a vital component in the 
study of the origin of resistance in disordered SC sys- 
tem, and will be used in a subsequent publication, to 
study the anomalous magnetoresistance observed in ex- 
periment 

We give a glimpse of such an analysis in the case of 
the superconductor-insulator transition in a disordered 
superconductor with increasing temperature. We take 
a 39% flUed model with weak disorder, set to = 0.2t, 
which displays a superconductor-insulator transition at a 
temperature ~ O.lAt. In Fig. 11 a) we show the varia- 
tion of conductance across the superconductor-insulator 
transition, and below it study the current distribution 
maps. In Fig. 11 'b) at T w 0.14Tc there are weak- 



disorder driven fluctuations in the SC order parameter, 
but an almost uniform supercurrent. The potential drops 
mainly in the contacts, and in the sample is equal to the 
average of the two leads with small random fluctuations. 
In Fig. [Tl|^d) at T « 2.3Tc the SC order parameter prac- 
tically vanishes, there is no supercurrent, and, due to 
the increasing resistance, only a small normal current 
flows through the sample. The potential, as expected 
for normal systems, decays linearly across the sample. 
At intermediate temperatures T Tc the current map 
Fig. [Tl^c) highlights the interplay of the normal and SC 
current. There is a rough correlation between regions of 
finite SC order parameter and supercurrent flow, on one 
hand, and zero SC order parameter and normal current, 
on the other. We point out three typical regions of the 
sample. Firstly, at (1) the order parameter is small and 
only normal current flows, whereas at (2) the order pa- 
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rameter is large and supercurrent flows. However, at (3) 
two SC regions are separated by a smaU normal region 
but are Josephson coupled and so a supercurrent flows 
through the zero SC order region. By examining the po- 
tential lines we see that the normal regions, for example 
(1), are acting as weak links whereas the potential drop 
over the superconducting regions is small. Thus the over- 
all resistance of the sample is dominated by such weak 
links. The current and potential maps allow us to see 
the superconductor-insulator transition developing, and 
we plan to investigate in details the relation of such a 
percolative picture to the Kosterlitz-Thoulcss transition, 
as was recently suggested [35]. 

V. DISCUSSION 

In this paper we have developed a new exact for- 
mula to calculate the current through a superconduc- 
tor connected to two non-interacting metallic leads with 
an imposed potential difference. The formula was im- 
plemented with a negative- 1/ Hubbard model which in- 
cluded both phase and amplitude fluctuations in the SC 
order parameter. A new Chebyshev expansion method 
allowed us to solve the model and calculate the current 
in 0{N'^^'^ M'^/'^) time, granting access to systems of un- 
precedented size. The formalism also enables the gener- 
ation of current and potential maps which show exactly 



where the super current and separately the normal cur- 
rent flows through the system. 

The formalism was exhaustively tested against a series 
of well-established results, demonstrating the accuracy 
of the procedure, its ability to capture various physical 
processes relevant to superconductivity in disordered sys- 
tems, and correctly model the presence of a magnetic field 
and finite temperature. These tests indicate that the for- 
malism and accompanying numerical solver can robustly 
calculate the current through a superconductor across a 
wide range of systems. In the future we plan to report 
on the application of the formalism to several outstand- 
ing questions, such as the magneto-resistance anomaly 
on crossing the superconductor- insulator transition |43) . 
the Little Parks effect in nano-scale cylinders [32], and 
dissipation-driven phase transitions in SC wires [44j . 
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Appendix A: Derivation of the Current Formula 

The formula for the current in the Bogoliubov basis set 

is 



■^=2|E/dKT^{t/L(e)r^-/R(e)r^] [u,(G>-G<)u*+v,(G>,-G<Jv;-av*(H>-H<)u*+au,(H>,-H<Jvj]} 

-fTr { - r«] [u*G<u* - v,G>^v: + au*H>v: - av,H<,u,] }) . (Al) 

We need to determine the Creen functions across the sample, which must be calculated in the presence of the leads. 
However, as the electrons in the metallic leads are non-interacting we can start from the bare electronic Green functions 
for the superconductor not coupled to the leads Gl^{m, n) = 5m,n/{^ — Cm + i<5) and G5j^(TO, n) = 5m,n/ (e + Cm + 
which have energy eigenstates Cm and (5 — >■ 0+ . We then write down Dyson's equation to self-consistently include the 
leads 

fQlA^Y2( ^lAKalpx'^P + Kaipx^p) ^^lai^gipxK - ^lalpxK) \ ( ^"A (A2) 

I 




Here 



IJL^ + \6) is the retarded 



Green function of the non-interacting electrons in the 
leads, with dispersion £p, and {up,Vp} are the matri- 
ces of the eigenstates multiplied by the lead plane wave 
states p at the tunneling barriers. To extract the re- 
tarded Green function and its anomalous counterpart 
from this matrix equation one has to perform a ma- 
trix inversion. The Dyson equation is for the retarded 
and advanced Green functions, whereas the current for- 
mula Eqn. (Al) is in terms of the lesser and greater 



Green functions. To transform these into the retarded 
and advanced Green functions we apply the identity 
G< = G< + G^Z^G< +G^Z<G- +G<Z^G- recursively to 
find G< = (1 -I- G^ZJ,)G<(1 + Z-G^) + G^5:<G^, where 
is the self energy. This recursion fixes the chemical 
potential of the superconductor by including tunneling to 
and from the leads. This will ensure that the net number 
of electrons is conserved, analogous to some extensions to 
the BTK formalism '271. However, as the final chemical 
potential must be independent of the chemical potential 
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of the uncoupled superconductor, the term containing 
must be identically zero leaving = G^^H'^G'^, and its 
greater Green function counterpart G> — GJ^Z^G^. We 



now extend this identity to include the anomalous Green 
function and recover 



H< 



CTH'^(up3< Vp - V. 



Up) H^(up5<p^u; + Vp5^^p^Vpj 



(A3) 



We can now take this, the analogous expression for the 
greater Green function, and their ano mal ous counter- 
parts, and substitute them into Eqn. (Al), which will 
yield Eqn. Q. 



Appendix B: Evaluation of the Monte Carlo 
Integrals 

In order to evaluate the correlation functions (e.g. 
Eq. 13), we need to sum over all possible spatial 



configurations of the auxiliary fields p and {A, A}, 
with each configuration carrying the weight P{p, A) ~ 
exp(— A])/Z. This distribution is sampled using 
the Metropolis algorithm [32|, which at each step pro- 
poses a new configuration of either the field p or A 
and calculates the resulting change in the total energy. 
If this change in the energy is negative the step is ac- 
cepted, whereas if positive it is accepted with probabil- 
ity exp{-f]{E[pncw] ~ E[poid])} and exp{-/3(i;[A„ow] - 
i?[Aoid])} respectively. Since the walk over p is one- 
dimensional we choose the step size |pncw ~ PomI to aim 
for 50% of the steps to be accepted, whereas the walk over 
{A, A} covers a two-dimensional space so we choose a 
step size | Anew — Aoidl so that 35.2% of the steps will be 
accepted 145'.. 

Central to the Monte Carlo method used to sample 
the partition function is the requirement to calculate the 
energy difference between two different configurations of 
the auxiliary fields, {poM, Aqm} and {p„ew,A„cw}- For 
a lattice with iV sites, to calculate the energy of each 
proposed configuration requires an effort of 0{N^), so 
an entire sweep over the N sites that make up the fields 
p and {A, A} requires a computational effort of 0{N'^). 
However, a recent method developed by Weifie |46j cal- 
culates just the difference between the energy of the con- 
figurations in a computationally efficient manner. For an 
update to the ith site a Chebyshev expansion with the 
< m < M coefficients containing {i\Tm{H /s)\i) must 
be calculated, where T™ is defined by the recursion re- 
lation T„(x) = 2xT„_i(x) - T™_2(x), To(x) = I, and 
Ti(x) = X. A typical expansion contained M — 1024 
terms. Previous authors have calculated this site- 
by-site through a succession of sparse matrix- vector mul- 
tiplications, each of cost 0{NM), so for an entire sweep 
over the order parameter the computational effort is 



0{N'^M). However, here we optimize the programme so 
that the entire sweep can be performed in 0{N^'^ M"^/^) 
time. Rather than follow a site-by-site approach cal- 
culated with sparse matrix- vector multiplications we in- 
stead calculate the matrix elements for the entire sweep 
simultaneously, which necessitates performing matrix- 
matrix multiplications. Provided the changes in the order 
parameters are small the local changes are independent 
of those of surrounding sites and we can then perform 
the entire sweep from this data set. Spherical averaging 
further reduces the influence of changes in the surround- 
ing order parameters. Central to the recursion relation 
for T,„ is the costly calculation of x", for \ < n < M . To 
evaluate this we divide the calculation of the M matrix 
products into three stages: 

1. The lowest order matrix products, up to x*"', are 
sparse. Therefore, for the elements 1 < n < k 
the matrix multiplications involve only sparse ma- 
trices, each of peak cost kN , and the total cost of 
calculating them is 0{k'^N). 

2. The second stage is to successively calculate every 
fcth matrix product. Each of these involves mul- 
tiplying the dense matrix x^^ by the matrix x'', 
for integer 1 < p < M/k, which costs ©(A^^.ss^ 
time HZj. With M/k of these products to calculate 
the total cost is 0{N^-^^M/k). 

3. The third stage is to construct the entire family of 
x" by interpolating between the matrices x^*^ found 
in the second stage. This is done by multiplying 
the dense matrices found in the second stage by 
the sparse matrices found in the flrst stage. Fur- 
thermore, as we need only the diagonal elements of 
the final matrix each separately costs 0{kN) and 
so the total cost is 0{kNM). 

Having now laid out the prescription of how to cal- 
culate the matrix elements, we now examine the total 
cost, 0{k^N + N'^-^^M/k + kNM). The choice k - 
^/W^^Hd will minimize the total cost to ©(iV^-^M^/^ + 
^1.46^4/3^^ and as typically N > M the cost is 
0(Ari-9Af2/3). This is a significant improvement over the 
cost O(N'^AI) of the Chebyshev expansion approach [46], 
which for the parameters employed in our simulations 
corresponds to a speedup by a factor of ~ 30. Now that 
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FIG. 12: (Color online) (a) The estimate of the current with 
number of Monte Carlo iterations, i, out of a total number 
/ = 1000. The primary y-axis shows the best estimate of 
the current (blue). The secondary y-axis shows the estimated 
standard deviation in this estimate (green) and idealized im- 
provement in the accuracy (red), (b) The distribution of 50 
separate current estimates at T = (red) and T = Q.2Tc 
(green) with best-fit Gaussian distributions, (c and d) The 
time r to perform a run on a 32 x 32 system renormalized by 
the time ro for a M = 512, A'^ = 1 system. In (c) the change 
with varying the system size A'^, where the blue line is for the 
standard 0{N'^) method of finding all of the energy eigenval- 
ues, the green is the 0{N^) standard Chebyshev expansion 
method g^, and the blue is the 0{N^-^) extended Cheby- 
shev approach. In (d) the two Chebyshev expansion method 
approaches are compared by varying the expansion order M. 
As in (c), the green line is the standard 0{A1) approach [46j . 
and the red line is the new ©(M^/^^) algorithm. 



the matrix elements behind the Chebyshev expansion 
have been found they are applied for the entire sweep. 

To verify the Monte Carlo procedure in Fig. [l2]^a) we 
first check the convergence of the estimate for the current 
and that its standard error falls as the root of the number 
of Monte Carlo iterations. In Fig. [l2jb) we compare the 
results of equilibrated Monte Carlo runs at zero temper- 
ature from a variety of initial configurations of the order 
parameter fields p and A. Evolution under the Metropo- 
lis algorithm drives these starting fields into different re- 
laxed configurations, which because the simulations are 
restricted here to T = are unable to be excited out 
to explore different configurations. These final configu- 
rations yield a variety of different current values, with 
standard deviation of ~ ±2.4% of the final total cur- 
rent. At finite temperature thermal excitations can drive 
the system to explore configurations around the ground 
state with a narrower standard deviation of ^ ±0.6%. 



Having verified the current statistics, in Fig. 12 'c and d) 
we show the results of some timing runs that highlight 
the improvement of the algorithm to 0{N^'^M^^^) time 
over the standard approach of calculating all the energy 
eigenvalues in 0{N^) time and the standard Chebyshev 
approach that runs in 0{N'^M) time. In particular, by 
varying the system size we observe that the method of 
calculating all the eigenvalues is more efficient for systems 
smaller than N ~ 10, but the new Chebyshev approach 
is superior for large systems. We took advantage of this 
development to study systems of unprecedented size. 

The Chebyshev expansion method just described rep- 
resents a zero order approximation. However, we can ex- 
tend this method further and calculate the lowest order 
change in the Chebyshev expansion following a shift in 
the configuration of the fields p and A by J. The resul- 
tant shift in the Chebyshev expansion of T^ is found using 
the recursion relationships t^ = |^Ti_i + |Hti_i — ti_2 
with to = and ti = 6/s. This allows the Chebyshev ex- 
pansion coefficients to be extrapolated over several con- 
figuration space sweeps, and the calculation time falls 
proportionally. Spherical averaging also reduces the in- 
fiuence of changes in the surrounding order parameters. 
In practice it was found that up to ten extrapolation 
steps could be performed, resulting in a code speed-up of 
a factor of ten. 

Though the Chebyshev approach can be used to di- 
rect the sampling of the system, to calculate expectation 
values, such as the current, it is necessary to diagonal- 
ize the system and determine the field configurations of 
its states. Formally this requires 0{N^) time. However, 
since the current is dominated by the quasiparticle states 
near to the Fermi surface we instead adopt the Implicitly 
Restarted Arnoldi Method [48] to calculate only those 
particular states. We are also helped by the sparsity of 
the matrix, which allows us to calculate K eigenstates 
in 0{KN) time. It is usually necessary to calculate a 
certain fraction of the energy states, so K (x N, and 
the total cost is 0(7V^). The eigenfunctions and energies 
can then be used to calculate the current for a specific 



17 



realization of p and A using the formalism described 
in Sec. It is then necessary to average over succes- 
sive realizations of p and A. However, the contribution 
from successive Monte Carlo calculations might be seri- 
ally correlated which would result in an underestimated 
value for the uncertainty in the predicted value of the 



current. To correct for this we calculated the correlation 
time through the truncated autocorrelation function [49j . 
Wc find a typical correlation time of approximately six 
Monte Carlo steps, which without autocorrelation correc- 
tions would correspond to an underestimate in the uncer- 
tainty of a factor of ^ 2.5. 
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